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ABSTRACT 

We demonstrate that the chaotic nature of iV-body systems can lead to macroscopic 
variations in the evolution of collisionlcss simulations containing rotationally sup- 
ported discs. The unavoidable stochasticity that afflicts all simulations generally causes 
mild differences between the evolution of similar models but, in order to illustrate that 
this is not always true, we present a case that shows extreme bimodal divergence. The 
divergent behaviour occurs in two different types of code and is independent of all 
numerical parameters. We identify and give explicit illustrations of several sources of 
stochasticity, and also show that macroscopic variations in the evolution can originate 
from differences at the round-off error level. We obtain somewhat more consistent 
results from simulations in which the halo is set up with great care compared with 
those started from more approximate equilibria, but we have been unable to eliminate 
diverging behaviour entirely because the main sources of stochasticity are intrinsic 
to the disc. We show that the divergence is only temporary and that halo friction is 
merely delayed, for a substantial time in some cases. We argue that the delays are 
unlikely to arise in real galaxies, and that our results do not affect dynamical friction 
constraints on halo density. Stochastic variations in the evolution are inevitable in 
all simulations of disc-halo systems, irrespective of how they were created, although 
their effect is generally far less extreme than we find here. The possibility of divergent 
behaviour complicates comparison of results from different workers. 
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1 INTRODUCTION 



Miller ( 19641 ) pointed out that all gravitational iV-body sys- 
tems are chaotic, in the sense that the trajectories of all 
particles in two systems that differ initially by a small shift 
in the starting position or velocity of even a single particle 
will diverge exponentially over time. Thus, two simulations 
started from the same initial conditions will follow identical 
evolutionary paths only if the arithmetic operations are per- 
formed with the same precision and in the same order, so 
that round off error is identical. These statements are true 
for every code, irrespective of the algorithm used for the 
computations, and no matter how many particles are em- 
ployed. In particular, a simulation can never be reproduced 
exactly when run with a different code. 

Microscopic chaos is unimportant for many applications 
because the different evolutionary paths of almost identi- 
cal simulations lead to similar macroscopic properties such 
as mass profiles, overall shape, etc., which therefore consti- 
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tute firm results. iBmnev fc Tremaind (|200Sl , her eafter BT08, 
344 ) make this argument and cite a test bv lFrenk et al\ 
1999) which indeed shows that many different codes yield 
similar key properties after following the collapse of a dark 
matter halo. In fact, results generally converge in tests that 
vary the n umerical grid, soft e ning, and/or number of par- 
ticles (e.g. iPower ef aZj |2003 : iDiemand et fflZ.ll2004l ). which 
they would not do if t here were a large element of stochas- 
ticity. ISellwoodl l|200Sl) also demonstrated exquisitely repro- 
ducible evolution of halo models that were perturbed by 
externally imposed bars, in sharp contrast to the results pre- 
sented here. 

Simulations with active discs of particles, on the other 
hand, are not so well behaved. ISellwood fc Debattistal (|2006l ) 
reported some minor differences, and one major, in a set of 
experiments using different numerical parameters but the 
same file of initial coordinates. We show here that simula- 
tions with discs can, at least for certain models, exhibit bi- 
modally divergent macroscopic results, even between cases 
that differ only at the round-off error level. The reason for 
this qualitative difference for discs is because collective insta- 
bilities and vigorous responses develop from particle noise. 
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Here we identify a number of distinct causes of stochastic 
behaviour in discs, and demonstrate explicitly how the evo- 
lution is affected. 

We show that the principal sources of divergent be- 
haviour are: (a) multiple in-plane global modes, (b) swing 
amplified noise, (c) bending instabilities, (d) suppression of 
dynamical friction, and (e) the truly chaotic nature of TV- 
body systems. We also show that the distribution of evolu- 
tionary paths taken in simulations of different realizations 
of the same model varies systematically with the care taken 
to set up the initial coordinates of halo particles. 

We deliberately choose to illustrate just how large the 
differences can be for one particular unstable equilibrium 
model. Stochasticity is present in all simulations and its ef- 
fects are always noticeable in those containing discs, but 
generally variations in the evolution show less scatter than 
in the case studied here. We show that the range of be- 
haviour is similar in two quite distinct TV-body codes and 
illustrate the sensitivity to differences at the round-off error 
level. We also show that increasing the number of particles 
does not reduce the spread of measured properties. 

Real galaxies are assembled and evolve in a compli- 
cated manner, and certainly do not pass through a well- 
constructed axisymmetric, equilibrium phase that is unsta- 
ble, although such a model is commonly used as a starting 
point of simulations. The objectives of experiments of this 
type are therefore (1) to determine whether plausible ax- 
isymmetric galaxy models are globally stable and (2) to de- 
velop an understanding of the dynamical evolution of mod- 
els that form bars and other non-axisymmetric structures. 
While we adopt a model of this type in this paper, its re- 
markable behaviour has implications for all simulations of 
disc-halo models, regardless of how they were created. 

The main part of the paper demonstrates the role of the 
five above-named sources of stochasticity in the evolution of 
disc models. We also explicitly show the effects of different 
particle selection techniques on the robustness of the be- 
haviour. Stochastic divergence has been reported elsewhere, 
but not recognize d as a n intrinsic aspect of these models; 
e.g., iKlypin et all {2008) attributed divergent evolution to 
inadequate numerical care, whereas stochasticity could be 
the cause. Appendix B reports extensive tests that confirm 
that the results we report here do not depend on any nu- 
merical parameters. 



2 SELECTION OF PARTICLES 

The selection of initial particle positions and velocities of an 
equilibrium model requires careful attention. Random selec- 
tion of even many millions of particles will lead to shot noise 
variations in both the density and velocity distributions of 
a model. Here we summarize the available techniques to se- 
lect initial coordinates of particles, with a focus on disc-halo 
models. These methods generally yield a set of particles that 
are not specific to any particular TV-body code. 



2.1 Selecting from a DF 

Jeans theorem requires that an equilibrium model should 
have a distribution function (DF) that is a function of the 
isolating integrals (BT08, p. 283). Thus the best way to 



realize an equilibrium set of particles for an initial model is 
to select from a DF, when one is available. 

While random selection of particles may be com- 
mon practice, it immediately discards a large part of 
this potential advantage. One w i dely used technique {e.g . 
Hollev-Bockelmann et al\ 120051: IWeinberg fc Katzl 120071 : 



Zhang fc Magorrianll200S ; iDubinski et aZ.ll2009h is to accept 



or reject candidate particles based on a comparison of a 
random variable with the value of the DF at the phase- 
space position of each particle, which introduces shot noise 
in the density of particles in integral space. The evolution 
of the simulation will be that of the selected DF, not the 
intended one, and different random realizations lead to sig- 
nificant variations in the m easured frequen cies of the insta- 
bilities in the linear regime l|SellwooJl983[ ) and substantial 
differences in the non-linear regime. It is therefore best to 
adopt a deterministic procedure for particle selection from 
a DF. 

A sche me to select par ticles smoothly in this way, 
first used in Sellwood 1 19831) a nd described more fully in 
ISellwood fc Ath anassoula ll986l) . is sum marized in the Ap- 
pendix of lDebattista fc Sellwood! |200oj). We divide integral, 
generally (E,L), space into n areas in such a way that 
J J FdEdL over each small area is exactly 1/nth of the 
integral over the total accessible ranges of 75 & L. Here 
F(E, L) is the differential distribution after integration over 
the other phase space variables (BT08, pp. 292, 299). Re- 
quiring that one particle lies within each area ensures that 
the selected set of particles is as close as possible to rep- 
resenting the desired particle density in integral space. We 
choose the precise position of a selected particle within each 
area quasi-randomly in order to ensure that the particles do 
not lie on an exact raster in integral space. We describe this 
scheme as deterministic selection from the DF, a term that 
ignores this minor random element. 

This scheme is readily adapted to select particles of un- 
equal masses if desired. To select particles having masses 
proportional to a weight function w(E,L), one simply 
weights the DF by w , which automatically adjusts the 
subdivision of ( E, Z/)-spa c e into areas of equal weighted DF, 
as described in lSellwoodl (|2008h . 

The phases of the particles around the orbit defined 
by these integrals can be selected at random. We have no 
evidence that the choice of radial phase, either for flat discs 
or for spheres, causes significant variations in the outcome 
and we discuss the choice of azimuthal phases in Section ^. 31 
below. 

iDebattista fc Sellwood! (|2000T ) describe the similar pro- 
cedure for 2-integral spheroidal models. 

2.2 When No Simple DF Is Available 

Comparatively few useful mass models have known DFs, 
and the realization of an equilibrium set of particles for 
a general model presents a significant challenge. Some au- 
thors (e.g. IShlosman fc Noguchi|[l993T ) have simply created 
a rough TV-body system, which they then evolve in the pres- 
ence of a frozen disc, thereby allowing the halo to relax to- 
war ds some ne a rby eq uilibrium. 

iHernquistl (|l993T l advocates solving the Jeans equa- 
tions for each component in the combined poten- 
tial of all mass components. His method is widely 
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used {e.g. Valenzuela fc Klypin 20031; lAthanassoulal 120031: 



lEl-Zant et al\\2Q04 ; iKlypin et aZ.I l2008h . but the resulting 
equilibrium is approximate. 

In general, it is better to derive an approximate DF for a 
spherical or spheroidal system. An isotropic DF for a spher- 
ical system can usually be obtained by Eddington inversion 
(BT08, p. 289), although it is important to verify that the 
function is positive for all energies (which it generally is, for 
reasonable mass models). 

Creating an equilibrium DF for a multi-component sys- 
tem presents a greater challenge, for which thre e effec - 
tive approaches hav e been develop ed. iRaha et al\ dl99lTi . 
iKuiiken fc Dubinskil (|l995l ) and iDebattista fc Sellwoodl 
1 20001 1 employ the method of iPrendergast fc Tomeij (|l970l ) 
to derive the mass distribution for a halo having some as- 
sumed DF that will be in equilibrium in the presence of one 
or more other mass components. Alternatively, one can use 
Eddington's inversion formula for the halo only in the poten- 
tial o f the combined disc and halo (|Hollev-Bockelmann et ali 
l2005l l. A third possibility, as here, is to start from a known 
spherical halo with a known DF and compress it by adding 
a disc and/or a bulge usin g Young's (1980) method (see 
ISellwood fc McGaughl 120051 ), and then to select particles 
from the compressed DF. Even though the last two methods 
use only the monopole term for the disc, all three methods 
yield a spheroidal system that is close to detailed equilib- 
rium everywhere. 

In general, it is more difficult to construct a good equi- 
librium for a disc component. The circular speed in the disc 
mid-plane as a function of radius is determined by the to- 
tal mass distr ibution and, commonly, one specifies Q(R) 
jToomrel H964I ) to determine the radial velocity spread at 
each radius. The Jeans equations in the epicycle approxi- 
mation (BT08, p. 326) generally yield a poor equilibrium 
except when the radial dispersion is a small fraction of the 
circular speed, and the asymmetric drift fo rmul a may have 
no solution near the centres of hot discs. IShul l|l969l ) de- 
scribes an approximate DF for a war m disc with a given 
radial velocity dispersion that we, and lKuiiken fc Du binski 
(1995), have found to be quite serviceable. Again in cases 
where the radial velocity dispersion stretches the validity of 
the epicycle approximation, radial gradients can lead to a 
disc surface density after integration over all velocities that 
differs slightly from that specified, as shown in Section [3. II 

The vertical structure of an isot hermal s t ellar sheet 
is given by th e formulae developed by ISpitzed l| 19421) and 
ICamml l|l950l ). and BT08 (p. 321) describe a generaliza- 
tion of the in-plane DF to include this feature, which they 
describe as the Schwarzschild DF. The Spitzer-Camm for- 
mulae assume full Newtonian gravity and no radial density 
or dispersion gradient. Force softening has an increasingly 
detrimental effect on the vertical balance as the ratio of disc 
thickness to softening length is reduced, we therefore pre- 
fer to construct a vertical equilibrium from the ID vertical 
Jeans equation in the actual force field of the softened disc 
potential, which leads to a better equilibrium. 

2.3 Quiet Starts 

The quiet start technique is a valuable addition to the set up 
process only when the model has a few vigorous, large-scale 
instabilities, such as arise in a cool, massive disc with a rota- 




Figure 1. The inner rotation curve of our standard model (solid). 
The separate contributions of the disc (dashed) and halo (dotted) 
are also shown. 



tion curve that rises approximately linearly from the centre. 
It is of little help when linear stability theory pr edicts the 
model to be responsive but ( almost) stable (e.g. ISellwoodl 
ll989l ; ISellwood fc Evaiisll200ll ). In these latter cases, collec- 
tive responses to residual noise grow more vigorously than 
any global modes, and the particle arrangement randomizes 
quickly. 

For a quiet start, one reproduces each selected master 
particle multiple times in a symmetrical arrangement, with 
image particles having the identical radius and velocity com- 
ponents in polar coordinates. We restrict the meaning of the 
phrase "quiet start" to this symmetrical arrangement of par- 
ticles - i.e. a, quiet start can be used no matter how the co- 
ordinates of the master particles are selected. Conversely, a 
"noisy start" means only that azimuthal coordinates are se- 
lected at random, again independent of how the master par- 
ticles are selected. The procedures for discs and spheroidal 
components differ slightly. 

For discs, we place image particles at the corners of an 
almost regular polygon in 2D, centred on the model centre. 
The polygon is not exactly regular because we nudge the 
particles away from exact n-fold symmetry by a random 
fraction of a small angle, typically 0.02°. When the disc has 
a finite thickness, the polygon must be duplicated with a 
second on the opposite side of the mid-plane for which both 
the vertical position z and velocity v z of every particle in 
each of the two polygons have opposite signs. 

When the force-determination method is based around 
an expansion in sectoral harmonics that is truncated at low 
order, m max , and the number of sides to the polygon n ^ 
2m max + 1, azimuthal forces in the initial model are much 
lower than would arise from particle shot noise - hence the 
label "quiet start". 

We have not tried quiet starts for other force methods, 
but they could still offer a significant advantage provided 
that the number of corners adopted for the polygon exceeds 
the azimuthal order of all the strong instabilities and non- 
axisymmetric responses (Section I5.3p by at least a factor 
two. 

We adopt a similar procedure for spheroidal compo- 
nents, except that we create image particles by rotating the 
initial position and velocity vectors using the usual rotation 
matrix for the adopted set of Euler angles (e.g. I Arfkenlll985l . 
p. 199). The set of Euler angles used creates an n-fold ro- 
tationally symmetric set of particles, which is also reflection 



© 2009 RAS, MNRAS OOO.IWlI 



4 J. A. Sellwood and V. P. Debattista 




Figure 2. Details of the approximate DF for the disc. Panels (a) and (b) show respectively the variation of / with radial velocity and 
azimuthal velocity at five different radii. Panel (c) shows the radial variations of the rms azimuthal speed (v^ solid) and radial speed (v r 
dashed), (d) compares the circular speed (dotted) with the mean v$ (solid) to illustrate the asymmetric drift. Panels (e) and (f) compare 
respectively the actual surface density and Q profiles (solid) with the desired profiles (dashed). The DF does not reproduce these curves 
perfectly, but the departures are minor. 



symmetric about the mid-plane, and has zero net momen- 
tum with a centre of mass at the model centre; each master 
particle is therefore inserted 2n times. It is reasonable to 
adopt n > 4. 



3 MODELS 

Here we describe all the various galaxy models we use in this 
paper. 



3.1 Standard Galaxy Model 

Our standard model is a composite disc-halo system with 
the rotation curve shown in Fig. [JJ The two mass compo- 
nents are an exponential disc and a compressed, strongly 
truncated, Hernquist halo. 

The initial surface density of the disc has the usual ex- 
ponential form 



E(fl) = 



M d 
2*< 



-R/Rd 



(1) 



where Md is the nominal disc mass. We truncate the disc 
at R — 5Rd, leaving an active disc mass of » 0.96M,j. The 
disc particles are set in orbital motion with a radial velocity 
spread so as to make Toomre's Q = 1.5. For most mod- 
els, we determine the approximate equilibrium velocities by 
solving the Jeans equations in the epicycle approximation 
as described in Section [272] 



In some cases we adopt Shu's approximate DF instead, 
and select disc particles deterministically from it. Properties 
of the DF and the radial variations of the low-order veloc- 
ity moments are shown in Fig. [5] While the radial velocity 
distributions are nicely Gaussian, the azimuthal velocity dis- 
tributions (J^Jj) are markedly skewed. This aspect, and the 
departures of the surface density and Q profiles from the de- 
sired values all decrease for models with less dominant discs 
or with lower values of Q. 

For fully 3D simulations, the density profile normal to 
the disc plane is Gaussian, with a constant scale height of 
0.05-Rd and appropriate vertical velocities in the numerically 
determined vertical force profile. 

We construct a halo in equilibrium with the disc in the 
following ma nner. We start fr om the initial density profile 
suggested bv lHernquistl (|l990h 



Po(r) = 



M h r s 



2irr(r a + r) 



(2) 



which has total mass and scale radius r s , with the 
isotropic distribution function (DF) also given by Hernquist. 
We strongly truncate this halo by eliminating all particles 
with enough energy to reach r > 2r s , causing the density 
to taper gently to zero at this radius, and an actual halo 
mass of w 0.25Mh- Since most of the discarded mass is at 
large radii, there is little change to the central attraction at 
r < 2r 3 and the model remains close to equilibrium. 

For our standard model, we choose r a = AQRd and set 
Mh = SOMd so that the halo mass is approximately 19 times 
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Table 1. Numerical parameters for our standard runs 



particle mass 

Figure 3. The frequency distribution of halo particle masses, in 
units of the disc particle mass. 



that of the disc. W e then employ the halo comp ression algo- 
rithm described bv lSellwood fc McGaughl (|2005l ) to compute 
a new, mildly anisotropic, DF for the compressed halo that 
results from including the above disc. The rotation curve, 
Fig. [1] shows that the disc dominates the central attraction 
over most of the inner part and the total rotation curve is 
approximately flat at large radii. 

We adopt a system of units such that G = Md — ad — 
1, where G is Newton's constant, Md is the mass of the 
untruncated disc, and ad is the length scale for the type of 
disc adopted. Therefore distances are in units of ad, masses 
are in units of Md, one dynamical time r = (aj/GAfd) 1 ' 2 , 
and velocities are in units of v = (GMd/ad) 1 ^ 2 = a d /r. One 
possible scaling to physical units is to choose the dynamical 
time to be 10 Myr and ad = 3 kpc, which implies Md = 
5.98 x 10 10 M s - The velocity unit v = 293 km s _1 , and the 
peak circular speed in Fig. [T]is approximately 235 km s~ . 

We also present results for two other disc-halo models 
for which we choose r s = 30R d and r s — 50Rd, i.e. that 
bracket our standard case. The more extended halo leads to 
a more dominant disc, while the disc is less dominant in the 
more concentrated halo. 

We select halo particles from the compressed DF using 
the smooth procedure summarized in Section [2.11 with the 
weight function for particle masses being w(L) = 0.5 + 20L, 
where L = \L\ is the total specific angular momentum. All 
disc particles have equal masses, but the masses of halo par- 
ticles range from 0.7 to 14.6 times the mass of the disc parti- 
cles. Fig. [3] shows the frequency distribution of halo particle 
masses. 

As a result of this careful procedure, both the disc and 
halo components are very close to equilibrium in the com- 
bined potential and the initial ratio of kinetic to the virial of 
Clausius (measured from the particles) is T/| W| = 0.498. At 
the same time, the phases of the particles in their carefully 
selected orbits are chosen at random, so that the model in- 
deed starts from the usual level of shot noise resulting from 
the random locations of the particles. 



3.2 Isochrone Disc 

We also present results using the isochrone disc with no halo. 
The potential (BT08, p. 65) has a simple form 
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Shortest time step 
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5 


5 



MR) 



while the surface density is 

E(i?) = 2^{ l0g ^ + (1 + a;2)1/2 ]-7 r 



2)1/2 



(3) 



(4) 



Here a is a length scale, and x — r/a; note E(0) = 
Md/(6iva 2 ). iKalnaisI l| 19761 ) describes a convenient family 
of DFs characterized by a parameter m,K\ we refer to each 
model as the isochrone/mx disc. He (jKal nais 1978]) also 
presents some preliminary results for th e normal modes, 
which were confirmed in simulatio ns dEarn fc Sellwoodl 
Il995l) . The local stability parameter (|Toomrdll964l ) for the 
isochrone/5 disc has an near constant value of Q ~ 1.6, and 
is Q ~ 1.2 for the isochrone/8 model. 



4 RESULTS 

We begin by showing just how much variation can occur. We 
first present the evolution of our standard disc/halo model 
whose rotation curve is shown in Fig. [T] Note that the disc 
equilibrium in these models is set by solving the Jeans equa- 
tions, while the halo particles are selected deterministically 
from a DF. Fig. U shows results from 16 separate runs with 
Sellwood's (2003) hybrid grid code using fixed numerical pa- 
rameters, given in Table [TJ but with different random seeds 
for the initial coordinates of the disc particles only. We plot 
the evolution of both the amplitude and pattern speed of 
the bar, measured as described in Appendix A. Even though 
the initial particles are selected from the same distributions, 
with different random seeds for the disc only, the amplitude 
evolution differs greatly from run to run and there is con- 
siderable spread in the evolution of the pattern speed. 

In order to demonstrate immediately that the scatter 
in Fig. U is not a numerical artefact of our grid code, Fig. [5] 
shows the results of a similar test with 5 runs using the 
tree code PKDGRAV l|StadeI 2001) using an opening angle 
8 = 0.7. PKDGRAV is a multi-stepping code, with time 
steps refined such that 5t = At/2' 1 < rj(efa) 1 ^ 2 , where e is 
the softening and a is the acceleration at a particle's current 
position. We use base time step At = 0.01 and rj = 0.2, 
which gives identical time steps for all particles. The results 
show a comparable spread in the evolution of both the am- 
plitude and pattern speeds. Results from the two codes with 
identical initial coordinates for all the particles do not com- 
pare in detail. For this problem, the tree code runs about 
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Figure 4. Evolution of the amplitude (left) and pattern speed 
(right) of the bar in 16 runs with different random seeds for the 
disc particle coordinates, run using Sellwood's (2003) hybrid code. 
The tiny differences in the initial models lead to a remarkably 
wide range of properties of the bar at late times. 




zoo 



Figure 5. Evolution of 5 runs with different random seeds for the 
disc particle coordinates, run using PKDGRAV with e = 0.05/?^- 



37 times more slowly than Sellwood's (2003) grid code; we 
therefore use it only for this cross check. 

The gross qualitative behaviour of all the models in 
Figs. [4] & [5] is similar at first. The bar forms at similar 
times with similar pattern speeds, though the initial peak 
amplitude varies by about ~ 25%. The evolution thereafter 
further diverges, notably with increasingly large differences 
in the bar amplitude. Steep declines in the bar amplitude 
in the interval 200 < t < 400 are gener ally associated with 
buckling events {e.g. iRaha et aZ.l ll991). but the timing of 
these events varies considerably. At late times in Fig. [4] the 
bar amplitude rises steadily in 9/16 simulations, although 
starting from different times in each case, while it stays low 
(over the time interval shown) in the remaining 7. 

It is more encouraging to note that the rate of decrease 
of the bar pattern speed does correlate with the bar ampli- 
tude; strong bars are more strongly braked by halo friction, 
as expected. Furthermore, continued amplitude growth of 
bars that are strongly b raked has been reported previously 
( e. a. lAthanassoulall2002l ) . 



4.1 Divergence at Late Times 

iDubinski et all (|2009h report a similar study of bar-unstable 
disc-halo models, which also reveal large amplitude differ- 
ences in the short term. However, they stress that the long- 




Figure 6. The inner rotation curves of models with (above) a 
slightly more domination halo and (below) a slightly more ex- 
tended halo. The line styles mean the same as in Fig. [TJ The 
behaviour of these models is shown in Figs. 171 fe 151 



term evolution of their simulations is reproducible, in con- 
trast to our finding. 

Fig. [7] shows that we confirm their conclusion for a dif- 
ferent model with a slightly more dominant halo; the evo- 
lution of both the bar amplitude and pattern speed shows 
much less scatter than is seen in Fig. [4] All cases show a 
steady rise in bar amplitude after the buckling event, al- 
though the curves for the different realizations during this 
stage of the evolution are offset in time, as also found by 
Dubinski et al. 

Fig. [8] shows results from a third model with a more 
dominant disc. The amplitude evolution in this model is 
again bi-modal, rising steadily at late times in half the cases, 
although not by as much as in our standard case (Fig. [4}. 
The rotation curves of both these models are shown in Fig. EI 

The late rise in bar amplitude occurs, if at all, only 
in models with live haloes and is associated with fric- 
tional braking. It is natural that frictional braking should 
be stronger when the halo is more dominant. In our stan- 
dard model (Fig. [4]|, and in the more dominant disc case 
(Fig. [8}, the large late-time differences arise because strong 
friction kicks in in some cases but not in all. We argue in 
Section [5. 5 1 that the reason for these differences is the exis- 
tence of adverse gradients in the halo DF, which can inhibit 
friction (|Sellwood fc Debattistall2006p . Whatever the cause, 
it is clear from these two sets of runs that onset of friction 
and steady bar growth at late times depends on compara- 
tively minor differences in the earlier evolution caused by 
the different random seeds. 

In order to quantify the scatt er, we compute the bi- 
weight estimate (|Beers et all Il990l ) of the mean and dis- 



© 2009 RAS, MNRAS 000,1111211 



Stochasticity in N-body Discs 7 




i — , — i — i — i — , — i — i 00 1 — i — i — i — i — , — i — i 00 1 — i — i — , — i — i — i — , — i 00 1 — i — i — i — i — i — l_ 

200 400 600 200 400 600 200 400 600 200 400 600 

time time shifted time shifted time 



Figure 7. Evolution of a set of models with a more dominant halo 
than those shown in Fig. [4] The initial rotation curve is shown in 
the upper panel of Fig. \E\ 
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Figure 10. Evolution of a set of 16 runs of our standard model 
that used a more careful disc set up procedure. 

0.4 I ' 1 ' 1 ' 1 ' 1 0.5 I ' 1 ' 1 ' i ' 








200 400 600 800 



200 400 600 BOO 



Figure 8. Evolution of a set of models with a less dominant halo 
than those shown in Fig. [4] The initial rotation curve is shown in 
the lower panel of Fig. \E\ 



Figure 11. Evolution of a set of 16 runs that used Hernquist's 
Jeans equation procedure to set up an approximate equilibrium 
for the halo particles. The bar amplitude grows at late times and 
the pattern decreases in all but three of these cases. 
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time time 

Figure 9. Comparison of the estimated means (solid lines) and 
±1<t scatter (dotted curves) in the three different haloes shown 
in Figs. H (red), [7] (blue), & [8] (green). 

persion of the measurements throughout all sets of exper- 
imentsQ Since bar growth is shifted slightly in time in the 
different runs shown in Figs. [4] [7] and [8] we apply a small 
time offset to the evolution of both quantities in order to 
ensure that the evolution coincides as the relative bar am- 
plitude grows through 0.1, before computing the mean and 

1 Their algorithm assumes the data to be unimodal with a few 
outliers, which is manifestly not the case in our data at late times. 



scatter from each set. Fig. [9] shows the time evolution of the 
means and scatter of the bar amplitude and pattern speed 
for all three haloes. It is clear that the stochastic spread is 
greatest for our standard halo (red lines), less for the less 
dominant halo (green lines) and least for the more dominant 
halo (blue lines). 

4.2 Particle selection 

Fig. llOl shows the consequence of selecting disc particles in a 
deterministic manner from an approximate DF as described 
in Sections 12.21 & 13.11 This procedure still has a random el- 
ement when choosing the precise values of E & L z within 
each sub-area, and the simulations have noisy starts because 
we randomly select the radial and azimuthal phases of the 
particles. The 16 different runs used different random seeds 
and are to be compared with those shown in Fig[4] for which 
disc particle velocities were selected from Gaussians whose 
widths were estimated from the Jeans equations in the epicy- 
cle approximation. There is no significant improvement, and 
in this case 6/16 runs have not slowed much by t = 800. 

The consequences of selecting halo particle velocities 
from Gaus sians whose widt hs are determined from the Jeans 
equations (|Hernquistll 19931 ) . are shown in Fig. 1111 With this 
more approximate halo equilibrium we see that all but 3/16 
bars grow and slow. The non-slowing fraction was 5/16 in 
a similar set of experiments (not shown) in which the halo 
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Figure 12. Evolution of the bar starting from 16 different selec- 
tions of particles from the same DF of the isochrone/5 disc. 



Table 2. Numerical parameters for our 2D simulations 





Isochrone 


Standard model 


Grid (N R ,N^) 


(180,256) 


(170,256) 


Sectoral harmonics 


^ m < 8 


^ m < 8 


Outer radius 


3.995a 


6.23i? d 


Softening rule 


Plummer 


Plummer 


Softening length e 


0.05a 


O.LRd 


Number of particles 


500 000 


various 


Equal masses 


yes 


yes 


Shortest time step 


0.025 


0.0125 


Time step zones 


1 


3 



particles were selected from the DF by the accept /reject 
method, instead of deterministically for Fig. [4] 

Thus we find a weak trend in these results with the 
quality of the different halo set-up procedures. The fraction 
of bars that do not experience strong friction rises to almost 
half when we use the most careful set-up procedure we have 
been able to devise for the halo, whereas use of the den- 
sity profile to choose radii and Jeans equations to set halo 
velocities results in a large majority (13/16) of bars that 
experience strong friction (Fig. Hip . This trend is also con- 
sistent with the weak dependence on halo particle number 
reported in Appendix B, where we find that the larger the 
halo particle number, the smaller the fraction of bars that 
slow. We also find a larger fraction of slowing bars when we 
use equal mass particles. These results hint that still larger 
calculations that are set up with extreme care may evolve 
in a consistent manner independent of the random seed, but 
we have been unable to demonstrate this. 



5 SOURCES OF STOCHASTICITY 

In this section, we describe and illustrate five sources of 
stochasticity, four of which contribute to the large scatter 
just described. 

5.1 A Reproducible Result 

We start from a simple unstable disc model for which the 
outcomes of simulations do not diverge with different ran- 
dom selections of initial particles. Fig. ll2l shows results from 
noisy start simulations in 2D of an isochrone/5 disc, in which 




shifted time shifted time 

Figure 13. The time evolution of the bar amplitude and pattern 
speed in a quiet start isochrone/8 disc in which Q ~ 1.2. Note the 
somewhat larger spread compared with that shown in Fig. 1121 
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Figure 14. Evolution of the bar in a noisy start isochrone/8 disc 
in which particles are drawn from the same DF as was used for 
Fig. [13] 



Q ~ 1.6; numerical parameters are given in Table [2] The dif- 
ferent curves come from separate simulations with different 
selections of particles from the same DF, using the "deter- 
ministic" procedure described in Section[5] The small scatter 
in the bar amplitude at late times can be further reduced 
by restricting disturbance forces to the m = 2 sectoral har- 
monic only. 

5.2 Multiple Modes 

Most unstable disc models support a large set of small- 
amplitude, unstaHejnodes_Juwinga wide range of growth 
rates (e.g. lToomrelll98ll ; iJalalil I2007T ). These linear modes, 
even those with the same angular periodicity, grow inde- 
pendently for as long as all disturbance amplitudes remain 
small. If the seed amplitudes of all modes are low, the first 
to saturate will be the most rapidly growing. In most unsta- 
ble discs, the fastest growing mode is generally the simplest, 
or fundamental, mode that is usually dubbed the bar mode. 
But if the growth rate of the bar mode does not exceed that 
of the next most vigorous mode by a large enough margin 
for some seed amplitudes, then both may have comparable 
amplitude when one saturates. The consequence of two or 
more modes reaching large amplitude at similar times but 
with random phases can lead to constructive or destruc- 
tive interference in the measured amplitudes as the "bar" 
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saturates. Non-linear effects then cause such differences to 
persist. 

We use the slightly cooler rriK = 8 isochrone disc to 
demonstrate this behaviour explicitly and, to avoid addi- 
tional complications, we restrict disturbance forces to those 
arising from the m = 2 sectoral harmonic only. Figs. [T3] & 
1141 illustrate the dependence of the outcome on the initial 
noise amplitude. The quiet start simulations in Fig. [13] are 
good enough that the growth rates of the two most rapidly 
growing m — 2 modes can be estimated by fitting to data 
from the extensive period of evolutio n before growth ends 
(e.g. ISellwood fc Athanassoura|[l986h . We find the growth 
rate of the second mode to be some 85% of that of the bar 
mode and that its amplitude (peak <5E/£) can be within a 
factor of a few of the dominant mode as the bar saturates. 
The consequence is a slight increase in the scatter of the 
later bar amplitudes in this case compared with the case for 
the hotter disc shown in Fig. 1121 

The mild scatter in Fig. [13] requires a quiet start, which 
decreases the seed amplitude of all non-axisymmetric dis- 
turbances that grow for ~ 100 time units before the rising 
amplitudes even become discernible in the figure. The much 
larger seed amplitudes when noisy starts are used do not 
allow the dominant mode to outgrow all others before sat- 
uration, with the consequences illustrated in Fig. 1141 The 
same sets of particles were used as for the results shown in 
Fig. 1131 but we placed the image particles at random az- 
imuths, instead of evenly. The period of rising amplitude is 
too short to allow more than very rough measurments of 
the growing modes, but it is clear that multiple unstable 
modes having comparable growth rates are seeded at large 
initial amplitudes by the shot noise. With such high seed 
amplitudes, there is not enough time for the most rapidly 
growing mode to outgrow the others, which therefore leads 
to very substantial variations in the final bar amplitudes. 
Note that this did not happen in the warmer disc (Fig. I12p . 
which also used a noisy start, since in that case all growth 
rates are lower, while the growth rate of the dominant bar 
mode exceeds that of all others by a larger margin. 

Notice also that not only is there greater scatter in both 
the bar amplitude and pattern speed in Fig. 1141 but both 
quantities scatter to lower values. We find indications that 
runs having lower pattern speed have the more dominant 
second mode. The fundamental bar mode, when it has time 
to outgrow the second mode, peaks at a greater amplitude 
and then relaxes back to lower value, as always happens in 
Fig. 1131 But when the second mode is competitive, the bar 
amplitude generally has a lower initial peak, and may even 
rise subsequently. 

5.3 Swing-amplified noise 

Our standard model is more complicated than the isolated 
isochrone disc. In particular, the inner rotation curve (Fig.[T| 
rises steeply where the halo density cusp dominates. Recall 
that a mode is a standing wave oscillation of the system, 
which can be neutral, growing, or decaying. The dominant 
linear global modes, known as cavity modes, in bar unstable 
discs are standing waves between the centre and corotation 
that must have a high en ough pattern speed to avoid any in- 
ner Lindblad resonances (|Toomrdll98ll ; lBmnev fc Tremainel 
120081 . pp. 508-518). The consequence of a steeply rising rota- 



tion curve is to make the maximum of the function £1 — k/2 
rise to high values near the centre, requiring any linear 
bisymmetric modes to have very high pattern speeds, small 
corotation radii, and to have very low growth rates (because 
the inner disc is not all that responsive). 

The outer disc, on the other hand, is highly responsive 
but has no cavity-type modes. We see evidence for weak 
edge-type mod es, which arise from a steep density gradient 
l|Toomrelll98l[ ) at the sharply truncated outer edge, but they 
are sufficiently far out and of low enough frequency to be 
decoupled from the bar forming process in the inner disc. 

Shot noise from the particles is vigorously amplified, 
but transient swing-amplified responses should be damped 
at the inner Lindblad resonance (ILR) of the disturbance 
l|Toomrelll98l ; lBinnev fc Tremainell2008l , p. 510), as long as 
the amplitude remains tiny. Large amplitude waves are not 
damped, however, and trap disc p articles near the ILR into 
a bar-like feature |Sellwoodlll989r i. 

Bar formation through amplified noise inevitably leads 
to a range of bar properties, but it is fortunate that the 
range turns out to be surprisingly narrow. To illustrate this, 
we study bar formation in our standard model in simplified 
simulations in which the motions of disc particles are con- 
fined to a plane, and the halo particles are replaced by a rigid 
mass component that simply provides the extra central at- 
traction to yield the same rotation curve as shown in Fig. [T] 
This approach has several advantages: the calculations are 
less expensive in computer time, but more importantly the 
dynamics is simpler because both bar buckling and halo fric- 
tion are eliminated, enabling us to isolate the bar formation 
process from these other complicating aspects of the overall 
evolution. 

Fig. [15] shows 4 sets of 16 runs each in which N is in- 
creased by a factor 10 from row to row, from N — 50K at 
the top, to N — 50M for the bottom row. The results from 
each run have been slightly shifted horizontally so that the 
amplitude passes through 0.1 at the same time (the mean 
for the 16 runs) as described above. The bar amplitude has 
a higher peak than in Figs. [4] & [5] in part, at least, because 
we use a different softening rule in 2D. The discrepant line 
in one of the pattern speed panels shows that the bar cannot 
always be identified in the early stages, but eventually it is 
in all cases. 

Fig. [T5] shows the evolution of the means and scatter 
in the four sets of experiments, and reveals that the main 
effects of increasing N are threefold: the formation of the 
bar is delayed because of lower seed noise, the mean peak 
bar amplitude increases and the scatter in the amplitude 
evolution rises with increasing particle number, at least to 
N — 5M. The pattern speeds are better behaved, with scat- 
ter decreasing as N rises. 

Because these calculations have less freedom, the am- 
plitude variation is much less than those shown in Fig. [4] 
which have the same numbers of disc particles as those in 
the second row of Fig. 1151 Nevertheless, the spread in the 
bar amplitudes after the initial rise remains quite high. The 
pattern speed does not decline as much because the rigid 
halo does not cause dynamical friction. 

Since amplified noise is intrinsically stochastic, the dom- 
inant transient responses in different random realizations of 
the disc must differ. The possible frequency range of the 
dominant pattern is broad, but not unbounded; the rota- 
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Figure 15. Evolution of the bar in four sets of 16 runs with dif- 
ferent random seeds for the disc particle coordinates. The number 
of particles rises by a factor of 10 from row to row, ranging from 
50K in the top row to 50M in the bottom row. 



tion curve and surface density profile, among other proper- 
ties, cause the responsiveness of the disc to vary with radius, 
and therefore the dominant responses have corotation radii 
in the region where the disc is most responsive. Thus the 
very first collective responses at low, but fixed, TV lead to 
initial bars having a range of strengths, i.e. sizes, with the 




200 400 " 200 400 



time time 

Figure 16. Summary plot showing the means and ±1ct scatter 
of the runs shown in Fig. 1151 



larger bars developing more slowly because the clock runs 
more slowly farther out in the disc. (The time delays have 
been removed from Fig. 1151 1 

The larger the number of particles, the longer it takes 
for the bar to form (Fig. ll6|) . Initial transient responses occur 
at roughly the same rate but, in experiments with larger 
TV, the lower initial amplitudes do not lead to immediate 
bar formation. Subsequent amplification events tend to be 
of greater amplitude, and to occur farther out in the disc. 
Thus we see that a lower level of shot noise favours large 
amplitude responses farther out in the disc that briefly lead 
to longer and stronger bars. 

The pleasant surprise is that after the initial transient 
episodes produce bars of different sizes and angular speeds, 
we observe (Fig. I16[) that subsequent evolution causes the 
range of bar strengths to narrow. Also most of the system- 
atic trends with particle number are erased in the subse- 
quent evolution, and neither the bar amplitude nor its pat- 
tern speed at later times exhibits more than a mild depen- 
dence on TV. It is fortunate that a degree of uniformity of 
the bar properties emerges after such tumultuously different 
evolution. But it is far from obvious why it should, especially 
since the model could have supported bars of wide range of 
sizes (e.g. Fig. [4} . 

The results shown in Fig. [15] are for models with rigid 
haloes in which the disc was created using the Jeans equa- 
tions (Section 12. 2|1 , Far from becoming better behaved, the 
scatter in the amplitude evolution increases as TV rises! We 
conducted a similar set of tests, also with rigid haloes, for 
which disc particles were selected deterministically from an 
approximate DF. The evolution of these more carefully set 
up models resulted in slightly improved behaviour: the bar 
formed somewhat more slowly, peaked at a little lower am- 
plitude for the same value of TV, and the scatter no longer 
varied systematically with TV. However, the final bar ampli- 
tude and pattern speeds were within the ranges shown in 
Fig.QU 

Unlike the results for the isochrone presented in Ap- 
pendix C, the more careful selection of particles yielded only 
a slight reduction in the spread in evolution. It is likely 
that this difference in behaviour of the two discs is due 
to the difference in bar forming mechanism; the instabil- 
ity of the isochrone disc is due to strongly unstable linear 
global modes, whereas as the bars in our standard model 
form through non-linear trapping of swing-amplified parti- 
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Figure 17. Comparison of the time evolution of two runs that 
differ only in the imposition of reflection symmetry about the 
midplane. The solid curves are for a model taken from Fig. [4] 
in which vertical forces are unrestricted while the dashed curves 
show the evolution of the same initial model when vertical forces 
from the disc are constrained to be symmetric about the mid- 
plane. 

cle noise that would be less affected by the quality of the 
equilibrium. 

Thus far we have discussed only bisymmetric instabili- 
ties, but other low-order instabilities may also be competi- 
tive. In fact, we find some evidence for lop-sidedness, which 
we describe in the next subsection. 

5.4 Bending modes 

The bars in most 3D simulations suffer from buckling in- 
stabilities that, when they saturate, thicken the bar in the 
vertic al direction (e.g. ICombes fe Sand ers 1981; Rah a et all 
1991). In many, but not all, cases the evolution of this bend- 
ing mode is quite violent and weakens the bar significantly, 
while the central density of the bar rises, as reported by 
Raha et al. The radial rearrangement of mass evidently lib- 
erates the energy needed to puff up the bar in the vertical 
direction. 

The time of saturation of the buckling mode depends 
on a variety of factors, such as the formation time of the 
bar, and the initial seed amplitude of the bending mode, 
the strength of the bar, etc. Several of these factors will in 
turn depend on the already stochastic formation of the bar. 
It is hardly surprising therefore, that this event occurs over 
a wide range of times and with a wide range of severity 
(Fig. [4]| , thereby compounding the overall level of stochas- 
ticity. 

The buckling mode can be inhibited by artificially 
imposing reflection symmetry about the mid-plane, which 
causes a substantial change to the evolution. Fig. [17] com- 
pares the evolution for one case; the dashed curves show that 
when buckling is inhibited, the bar continues to grow in am- 
plitude, while slowing, for a long period. On the other hand, 
the amplitude drops quite abruptly when the bar buckles 
(solid curves) and the subsequent amplitude and pattern 
speed hold approximately steady. 

Not all the bars in the runs illustrated in Fig. [4] experi- 
ence a violent buckling event. In some cases the bar ampli- 
tude does not decrease after the initial peak, while in others 
the amplitude drop is more gradual. 

Fig. [18] shows the effect of suppressing the m = 1 sec- 
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Figure 18. Evolution of a set of 16 runs that differ from those 
shown in Fig. [4]only in the suppression of lop-sidedness about the 
z- axis. 

toral harmonic about the z-axis for both the disc and halo 
particles. This has the effect of preventing the centres of ei- 
ther component from leaving the z-axis. (Suppressing the 
/ = 1 component of the halo force calculation would nail the 
centre of that component to the origin, which would pre- 
vent the halo from responding properly to a buckling mode.) 
With lop-sidedness inhibited in this way, all bars buckle, and 
all but one do so violently with a large decrease in ampli- 
tude. This difference in buckling behaviour from that shown 
for the same initial models in Fig. [4] indicates that buckling 
is strongly influenced by mild lop-sidedness, which has not 
been reported elsewhere, as far as we are aware. We could 
not find any evidence for lop-sided instabilities in the runs 
shown in Fig. [5] and the distance between the centroids of 
the halo and disc particles was < 0.002^. As it seems un- 
likely that such small offsets could have such a large effect 
on the saturation of the buckling mode, we think it possible 
that an anti-symmetric mode competes. Investigation of this 
possibility here would be too great a digression. 

Despite the violence of most buckling events, most bars 
in these restricted simulations continue to slow after the 
buckling event and amplitude growth resumes. The four ex- 
ceptions are bars that remained strong right after their for- 
mation and did not slow much either before or after the 
buckling event. 

Results reported in Appendix B show that the buckling 
beh aviour is also s omew hat sensitive to particle softening. 

iKlypin et al. I (|2008h report that the violence of the 
buckling event also depends on the initial thickness of the 
disc. This is as expected, since iMerritt fc Sellwoodl (|l994T ) 
showed that buckling is a consequence of a collective in- 
stability that arises in systems in which the velocity dis- 
tribution becomes too anisotropic, and thickening the disc 
reduces the flattening of the velocity ellipsoid. However, in a 
separate test with a set of runs with twice the disc thickness 
(not shown), we still find a similar degree of scatter in the 
late evolution. 

5.5 Incidence of Dynamical Friction 

Fig. [19] shows that the divergent late-time evolution of the 
runs shown in Fig. [4] is due to differences in the incidence 
of dynamical friction. The lines are coloured blue when the 
torque acting on the halo dL z /dt < 5 x 10~ 5 GM 2 / Rd, and 
are red otherwise. 
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Figure 19. The results shown in Fig. [4] but with the curves 
colored blue when the torque on the halo is low and red otherwise. 
The calculations were continued for models that had not slowed 
by t = 1000 and were stopped either at t = 3000 or soon after 
friction kicked in, which happened in all but two cases. 

The absence of bar friction may have a variety of 
causes: (a) low halo density, (b) a weak bar, and (c) 
metastability caused by local adverse gradients in the den- 
sity of halo particles as fun ction of angular momentum 
ijSellwood fc Debattistal l2006h . The halo density is just 
about the same in all cases, but the bar strength varies 
widely and it is clear that the weaker bars experience lit- 
tle friction. 

The third possibility is indicated by the evidence in 
Fig. 1191 since friction eventually resumes, sometimes after 
a very long period during which the bar amplitude does 
not increa se; the metastable state does not last indefinitely. 
We argue (|Sellwood fc Debattista|[2 006) that the metastable 
state has a finite lifetime because weak friction at minor 
resonances gradually slows the bar until the more impor- 
tant resonances move out of the region of adverse gradients, 
allowing strong friction to resume. 

Metastability could be caused by the buckling event, 
since bars that are weakened substantially by a buckling 
event, such as the case picked out in Fig. 1171 generally do 
not experience much friction at late times, and their ampli- 
tudes stay low. The upward rise in the bar pattern speed 
at the time of bucklin g is shown clearly by the soli d curve 
in Fig. [171 which we ijsellwood fc Debattistal 120061 ) found 
to be a likely cause of metastability. It is reasonable that 
the concentration of mass to the centre as the bar buckles 
should cause an upward fluctuation in the bar pattern speed 
(because the orbit periods must vary inversely as the square 
root of the mean interior density). However, buckling does 
not always lead to a cessation of friction; for example, many 
of the bars in the 16 runs with a more dense halo (Fig. [7} 
clearly buckled, but friction continued in all cases. 

5.6 True chaos 

Here we show that Miller's (1964) instability can lead to 
macroscopic differences in discs. Where initial evolution is 
largely determined by swing-amplification of the spectrum of 
particle noise laid down by the random coordinates of par- 
ticles, models that differ by tiny amounts quickly diverge 
because the subsequent spiral events depend on the details 
of evolution of previous events. This phenomenon causes the 
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Figure 21. The upper panels compare the evolution of four cases 
that started from the identical file of particle coordinates, with 
all numerical parameters held fixed, except that solid lines are for 
calculations in single precision, dotted lines are for the identical 
calculations in double precision. As for Fig. 1201 the order of the 
particles was reversed in one of each pair. The lower panel shows 
the time evolution of the quantity d defined in eq. (0 for both 
pairs of runs. 

micro-chaos in iV-body systems to lead to macroscopic dif- 
ferences in discs. 

Fig. [20] compares the amplitude evolution of each case 
shown in Fig. U (solid lines) with another run of the same 
case with the order of the particles reversed (dashed) . Thus 
the initial phase space coordinates of all particles were iden- 
tical and are evolved with the same code on identical pro- 
cessors. Each pair of simulations differ only in the order in 
which arithmetic operations are performed, which changes 
the initial accelerations at the round off error level only, yet 
the amplitudes at late times generally differ visibly, and in 
some cases, e.g. 10 & 15, the evolution differs qualitatively. 

So far, every calculation with grid codes that we have 
reported here was conducted using single precision arith- 
metic for most operations. We have checked that increased 
precision has no effect on the range of behaviour shown in 
Fig. [4] and results differ only slightly, as we now show for 
one case. 

Fig. [21] shows that the system remains chaotic when 
we repeat the calculations using double precision arithmetic 
(dotted lines). The higher precision calculations begin to 
diverge visibly at about the same times as in the single pre- 
cision cases, and the subsequent differences are comparable. 
In order to monitor the divergence in these cases, we com- 
pute the value over time of the difference 

d = [&(A2, a - A 2 , b f + S(A 2 , a - A 2 , b f] 1/2 (5) 
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Figure 20. Comparison of the amplitude evolution of the models shown in Fig. [4] (solid lines) with the same sets of particles processed 
in reverse order (dashed lines). The evolution of these two sets of identical runs is measurably different in all cases, and qualitatively 
different in some, especially cases 10 & 15. The dotted lines in the first 5 panels show the evoltions using PKDGRAV for the same files 
of initial particles. 



between the bar coefficients (eq. I A3f) in these pairs of exper- 
iments (a & b) in which the order of the particles was re- 
versed. The solid (dotted) line in the lower panel of Fig. 1211 
shows the result for the single (double) precision pair. By 
t ~ 300 the models differ quite visibly in the amplitude and 
phase of the bar, which accounts for the fact that d asymp- 
totes to a lasting value where the phases of the two bars 
differ. 

The difference, d, in double precision grows quasi- 
exponentially over time at first, which is symptomatic of 
chaos, with a Lyapunov (e-folding) time of ~ 4.75 dynami- 
cal times, i.e., less than 25% of the orbit period (~ 20 dy- 
namical times) at R — 2.5Rd- Using this estimate of the 
Lyapunov time, the difference in the double precision case 
should equal the initial difference in the single precision case 
after m 93 dynamical times, and the early evolution of d in 
the lower precision case is roughly similar to that in the 
double precision case with a time offset of this magnitude. 
Even though there is a much smaller initial difference be- 
tween the two double-precision models, the seed amplitude 
of the instabilities is set by the shot noise, which is the same 
in all 4 runs. Thus the non-axisymmetric structures are al- 
most fully developed in the double precision models by the 
time the dotted curve reaches the level of the start of the 
solid line; therefore one cannot expect the curves to overlay 
perfectly. 

It is curious that the difference in the double precision 
case "catches up" with that in the single precision case. The 
shoulder in log 10 d that appears in both precisions at about 
t = 300 seems to be responsible for this convergence, which 
occurs both at such a large value of d as to be well past where 



exponential divergence could be expected to hold, and at a 
time when the bar in all four runs is fully developed. 

A perfect collisionless particle system should be exactly 
time reversible; that is, if the velocities of all the particles 
were reversed at some instant, the system should retrace its 
evolution. Fig.[22]shows that reversed simulations do retrace 
their evolution for a short while, between 60 and 80 dynami- 
cal times, after which the evolution of the reversed model vis- 
ibly departs from the corresponding reflection of the forward 
evolution. This period of successful reversibility is consistent 
with our Lyapunov divergence estimate: 15 Lyapunov times 
(= 71.25 dynamical times) corresponds to a divergence of 
~ 10 6 ' 5 , which is sufficient to alter almost every significant 
digit in these single precision calculations and lead to re- 
versed evolution that becomes largely independent of that 
in the forward direction. Further analysis of these simula- 
tions revealed that the first signs of irreversibility appeared 
as differences in the leading spiral Fourier components, sug- 
gesting that vigorous swing-amplification of particle noise is 
primarily responsible for the short Lyapunov time. 

We conclude from these tests that the A-body system 
we are trying to simulate is indeed chaotic. Further, the 
effects of chaos are not significantly worsened by round-off 
error in single precision; we have also verified that the full di- 
vergence of the results in Fig. [4] persists in double precision. 
In fact, the first author has frequently checked, and always 
confirmed, that no advantage results from use of higher pre- 
cision arithmetic when computing the evolution of collision- 
less A^-body systems. This conclusion is in sharp contrast 
with the requirements for collisional systems (e.g. lAarsethl 
l200St) . 
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Figure 22. The magenta line shows the unsmoothed bar ampli- 
tude evolution of one model run to t = 200. The other lines show 
the continued evolution of the same model with the velocities of 
all the particles reversed at t = 50 (cyan), t = 100 (blue), t = 150 
(green), and t = 200 (red). In all four cases, the evolution imme- 
diately after the reversal faithfully retraces the forward evolution 
for a period less than 100 time units. After this time the evolution 
departs noticeably from a reflection of the line about the reversed 
moment. 

In none of the simulations with grid codes reported in 
this paper did we distribute the computation over multiple 
parallel processors, even though the code has been well op- 
timized for parallel use. We adopted this strategy in order 
to avoid the additional randomness that is inevitable when 
results from multiple processors are combined in an unpre- 
dictable order. 

The dotted curves in the first five panels of Fig. I20l show 
the result using the tree code PKDGRAV for the same ini- 
tial coordinates in each case, which are reproduced from 
Fig. [5] Although the ranges and distributions of measured 
bar properties shown in Figs. [4]&[5]are similar, the results 
do not compare in detail, as noted above. Results from the 
two different codes diverge strongly in all but one case, re- 
inforcing the conclusion of intrinsic stochasticity. Which of 
the two possible evolutionary paths is taken in the evolution 
is affected no more, and no less, by code differences than by 
choices of the random seed. 



6 DISCUSSION 

6.1 Is there a right result? 

One of the most troubling aspects of the diverging evolution 
in Figs. |3] & O is that one cannot decide which of the two 
patterns of behaviour is "correct," or indeed whether there 
could be a unique evolutionary path with a perfect code and 
infinite numbers of particles. 

Since these models have high density centres (Fig. [TJl, 
linear stability analysis would most likely reveal that all 
global modes, with the possible exception of edge modes 
(Toomrej[T98lh . have very low growth rates, and therefore 
the disc ought to be stable and not form a bar. If this is 



indeed what linear theory would predict, then the "right re- 
sult" with a perfect code and infinite numbers of particles 
would be a stable model that does not form a bar. This 
outcome never occurred in the > 400 simulations we report 
here, even in cases with one hundred times our standard 
number of disc particles (Fig. I15|l . 

The level of shot noise in a simulation with > 1 mil- 
lion particles is clearly ~ 100 times higher than would be 
present in a real galaxy if the ~ 10 10 stars were randomly 
distributed. But the mass in real galaxy discs is clumpier be- 
cause of the existence of star clusters and giant gas clouds, 
which raises the amplitude of random potential fluctuations 
- although the density fluctuation spectrum may not be the 
same as that of shot noise in the simulations. Nevertheless, 
it seems most unlikely that a real galaxy closely resembling 
the model used in our simulations could avoid being barred. 

6.2 Dynamical Friction 

The greatest source of divergence is the bimodal nature of 
dynamical friction, which is avoided for a long time in some 
cases, but kicks in immediately in others, causing the bar 
to slow and increase in strength by a substantial factor. It 
is likely that friction is avoided because the needed gra- 
dient in the halo DF as a function of angular momentum 
has been flatt ened by the earlier evolution of the model, 
as reported bv lSellwood fc Debattista! (|2006l ). The fact that 
this happens her e more frequently than we f ound with the 
model created bv lValenzuela fc Klypinl (|2003l ) may have two 
causes: their model had both a less dominant disc and an 
initial halo with significant departures from equilibrium. 

In Section 14.21 we reported a weak trend towards a 
larger fraction of non-slowing bars as we took greater care 
over the initial selection of particles; further, the largest frac- 
tion (10/16) occurs in the test with four times the number 
of halo particles reported in Appendix B. This weak trend 
suggests that the metastable state is reached more readily 
as the quality of the simulation is improved . 

However. [Sellwood fc Debattis"taf (|2006h found that the 
metastable state, in which the bar did not slow, was not in- 
definite and friction eventually resumed, as we also find here 
(Fig. I19[l . Furthermore, they found the metastable state to 
be fragile, and friction would resume soon after a tiny per- 
turbation, such as the distant passage of a small satellite 
galaxy. Thus, even though the metastable state is reached 
more frequently in higher quality calculations, it is unlikely 
it could be sustained in real galaxies. We conclude there- 
fore that the strongly braked and growing bar is the most 
"realistic" outcome from these simulations. 

6.3 Introducing a seed disturbance 

Irlollev-Bockelmann et al\ (|2005t ) attempted to make the out- 
come more predictable by seeding the bar instability by an 
externally applied transient squeeze. We argue here that this 
approach is not the panacea it may seem. 

In the case of discs having well-defined global instabili- 
ties, noisy starts already seed t he dominant u nstable modes 
at high amplitude (Section l5. 21 ISellwoodll 19831) . If a seed dis- 
turbance is to prevail, it must be imposed at such high am- 
plitude as to be practically non-linear at the outset. Further- 
more, the objective must be to favour the dominant mode 
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over the others, which cannot be achieved by a simple per- 
turbation. Instead, one must impose both the detailed radial 
shape and perturbed velocities of the mode, which are gen- 
erally not known. A more generic disturbance, such as a 
"squeeze" will simply raise the amplitude of all the modes 
and transients, giving less time for the do minant mode to 
outgrow the others. Quiet starts (Section [2.3llSellwoodlll983l ; 
ISellwood fc Athanassouiall 19861 ). however, have the effect of 
reducing the initial amplitudes of all non-axisymmetric dis- 
turbances to such an extent that there is ample time for the 
most rapidly growing mode to prevail. Thus the outcome of 
a quiet start experiment is tolerably reproducible without 
the need to apply an additional seed (Fig. 1 13 p . 

The situation is far more difficult in the case, as in the 
present study, where the disc has no prevailing global insta- 
bilities, since the evolution of a simulation is dominated by 
swing-amplified shot noise. Quiet starts are all but useless in 
these circumstances also, since they break up rapidly as the 
tiny seed noise is swing amplified, with similar outcomes, 
only slightly delayed, to those from noisy starts. Cranking 
up the particle number does not reduce variations in the bar 
amplitude at later times (they actually increased in Fig. ll5[) . 
but does delay bar formation. Because of this, perhaps a 
suitable seed disturbance in a very large iV disc may prevail 
over the amplified shot noise and lead to a more reproducible 
outcome. We have not explored this idea here and leave it 
for a future study. 



7 CONCLUSIONS 

We have shown that simulations over a fixed evolutionary 
period of a simple disc-halo galaxy model can vary widely 
between cases that differ only in the random seed used to 
generate the particles, even though they are drawn from 
identical distributions. Fig. [3] shows that the late-time am- 
plitude of the bar can differ by a factor of three or more while 
the stronger bars may have half the pattern speed of the 
weaker ones. Fig. \T§\ shows that the largest differences are 
only temporary, however. We have deliberately focused our 
study on a case which displays this extreme bad behaviour. 
Stochastic variations are inevitable, but evolution is gener- 
ally less divergent; e.g., when the halo has both higher and 
lower density {e.g. Fig. |5J. 

We have shown that the divergent outcomes do not re- 
sult from a numerical artefact, since they are independent of 
numerical parameters (Appendix B). Also, similar behaviour 
occurs with a code of a totally different type (PKDGRAV, 
see Fig. [5j. Instead, this extreme stochasticity results from 
a number of physical causes that we have identified and il- 
lustrated. The most important for our model are: swing- 
amplified particle noise, the variations in the incidence and 
severity of buckling, and the incidence of dynamical friction. 
We have separately shown (Fig. I14|l that other disc models 
having a well-defined spectrum of global modes can have a 
range of outcomes because of the coexistence of competing 
instabilities. 

The calculations in Fig. [4] are of models that were set 
up with considerable care so as to be as close as possible 
to equilibrium. An additional level of unpredictability can 
result from less careful set-up procedures, as illustrated in 
Appendix C. 



We have been aware for many years that simulations 
including disc components can be reproduced exactly only 
if the arithmetic operations are performed in the same order 
to the same precision, and that differences at the round-off 
error level can lead to visibly different evolution. However, 
we have been surprised by the strongly divergent behaviour 
of the particular model studied here. The pairs of divergent 
results in Fig.[20]are the stellar dynamical equivalents of the 
possible macroscopic atmospheric consequences of Lorenz's 
butterfly flapping its wings. Because the system is chaotic, 
improved precision arithmetic is of no help in reducing the 
scatter in the outcomes. 

The divergence in different realizations of our standard 
case arises from a temporary delay in the incidence of dy- 
namical friction, which is determined by minor details of 
the early evolution. Strong friction causes the bar to both 
slow and grow; in some cases this occurs right after bar 
formation, but in others the bar rotates steadily at an al- 
most constant amplitude for a protracted period. Friction 
is avoided when the earlier evolution causes an inflexion in 
the angular momentum dens ity gradient of the halo. We 
(|Sellwood fc Debatti sta 2006) previously described this as 
a metastable state because it did not last indefinitely even 
when the evolution was unperturbed, and we also showed 
that mild perturbations could cause friction to resume. We 
find that the fraction of initially non-slowing bars increases 
as greater care is taken over the initial set up because the 
smaller fluctuations in such models are less likely to nudge 
the model out of the metastable state. 

We argue in Section [S] that the most realistic outcome 
of these experiments is the slowing and growing bar, despite 
the fact that we find the delayed friction result increasingly 
often as we improve the quality of the initial set-up and 
of the simulation. Since most real galaxies are likely to be 
subjected to frequent mild perturbations, we conclude that 
slowing and growing bars are in fact the more realistic out- 
come. 

Since the possible evolution of the simulation is not 
unique, multiple experiments of essentially the same model 
are needed in order to demonstrate that the behaviour is 
robust. Furthermore, the failure of an experiment by one 
group to reproduce the results of a similar experiment by 
another may not be the result of errors or artefacts in ei- 
ther or both codes, but rather a reflection of a fundamental 
stochasticity of the system under study. 

iKlypin et all (12008T ) report a similar, but less extensive, 
comparison between two tree codes and an adaptive mesh 
method, and conclude that all the codes produce "nearly 
the same" results in simulations performed with sufficient 
numerical care. However, inspection of the comparatively 
short evolution shown in their Fig. 8 reveals slowly diverg- 
ing outcomes, even between two simulations run with tree 
codes. They also report (their Fig. 1) a strongly divergent 
result when the time step was varied; the sharp decrease in 
bar strength in this one case was clearly a consequence of a 
more violent buckling event than in their comparison cases. 
Such a difference could have easily arisen from stochastic 
variations of the kind discussed here, and the conclusion 
that the shorter time step is re quired no longer follow s. We 
show here (Appendix B), as do lDubinski et all (j2009T ). that 
results are robust to wide variations in time step. Clearly 
when stochasticity can lead to sharply divergent results, pa- 
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rameter tests that throw up surprises are conclusive only 
after ensembles of particle realizations have been simulated. 
This must also be a requirement for meaningful comparisons 
between codes or workers. 

Since the principal sources of stochasticity are con- 
nected to disc dynamics, they are unrelated to the h a lo par - 
ticle number question r aised by IWeinberg fc Kat j (|2007h . 
Not only has ISellwoodl | 2008) already shown that friction 
can be captured adequately with moderate particle num- 
bers, but we have found here that the expected bar friction 
arises more readily in haloes with fewer or equal mass halo 
particles, or in haloes that are not set up with great care - 
which is not the expected behaviour were particle scattering 
dominant. Instead, small departures from equilibrium can 
upset the delicat e metastable state in which b ars can rotate 
without friction (|Sellwood fe Debattistall2006l ). 

It should be noted that bars that slow through dy- 
namical friction also grow in length, as reported earlier by 
lAthanassoulal l|2002h . Nevertheless, for these models the ra- 
tio of corotation radius to bar semi-axis 1Z > 1.4, as expecte d 
for a moderate-density halo (|Debattista fc Sellwood! 12000). 
Those bars that avoid frictio n for a long period, however , 
have 1Z < 1.4, as also found bv lValenzuela fc Klypinl l|2003l) . 
but this met astable state is fragile and un likely to arise in 
real galaxies l|Sellwood fc Debattista|[200^ ). 

Since all iV-body simulations are intrinsically chaotic, 
they can be reproduced exactly only if the same arithmetic 
operations are performed in the same order with the same 
precision, as noted in the introduction, and borne out in 
Fig. [20] These requirements dictate the use of the same 
code, compiler, operating system, and hardware. Further, 
if the calculation is stopped and then resumed, it is impor- 
tant to save sufficient information so that the acceleration 
used to advance each particle at the next step is identi- 
cal, to machine precision, to that it would have been had 
the calculation not been interrupted. This can be arranged 
without too much difficulty, if the calculation is run on a 
single processor. However, simulations that distribute work 
over parallel processors in computer clusters would be ex- 
actly reproducible only if care is taken to ensure that the 
work is distributed and the results are combined in a fully 
predictable manner. 

Provided the divergence is slight, exact reproducibility 
is of little scientific interest, although such a capability is 
useful to the practitioner. But when, as described here, the 
model under test can have strongly divergent behaviour that 
arises from differences that begin at the round off level with 
the same code on the same machine, comparison of results 
between different codes and on different platforms becomes 
much less likely to produce agreement, even when the simu- 
lations share the same file of initial coordinates. It is ironic 
that the model used here was in fact that selected as a test 
case for code comparison; fortunately, the authors discov- 
ered its unsuitability in time! 
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APPENDIX A: CODES AND SOFTENING 
RULES 

Al Force Determination Methods 

The accelerations to be applied to particles in an iV-body 
simulation can be determined in many different ways that 
fall into two broad classes. Direct pair wise summation, usu- 
ally with a tree algorithm to improve efficiency, and methods 
that solve for the gravitational field over a volume. Three 
common methods in the latter category are: (1) solving a 
finite difference approximation to the Poisson equation on 
a grid, (2) convolution between the source distribution and 
a Green function on a grid, and (3) expansion of the field 
in multipoles, with either a basis set to represent the radial 
part or a grid on which the contributions of interior and ex- 
terior masses are tabulated. Grid and field methods are far 
more efficient than tree codes, albeit at the cost of ease of 
use and versatility. 

All grid methods assign masses to a spatial raster of 
points and tabulate the field at the same points. Sensible 
interpolation schemes to treat particles between grid points 
lead to forces between particles that decrease smoothly at 
separations below one grid space, reaching zero for coinci- 
dent particles. 

Finite difference methods solve an approximation to the 
Poisson equation directly, yielding a potential arising from 
the mass distribution. Acceleration components, which have 
to be estimated from a finite difference approximation to the 
gradient operator, lead to forces that approximate the full 
Newtonian value at distances of greater than a few mesh 
spaces, but whic h are significantly weake r at short range 
(e.g. appendix of lSellwood fc Merritd[l994h . 

Convolution methods, on the other hand, can be used 
to compute the acceleration components directly, without 
the need to difference a potential. The Green function is the 
force field of a unit mass, which requires a separate convo- 
lution for each coordinate direction. However, the force law 
needs to be softened at short range both to prevent accel- 
eration components from varying so steeply across a grid 
cell that simple interpolation rules become inadequate and 
also to limit the maximum possible acceleration, particularly 
where grid cells become very small near the centres of polar 
grids. 



A2 Softening Rules 

Since any arbitrary softening rule can be adopted in convo- 
lution methods, physical considerations can be used to select 
the optimum rule for a particular application. The Plummer 
softening rule for a unit mass uses the density profile and 
potential 

p(*) = 4^( l +* 2 r 5/2 ; = -7(i+z 2 r 1/2 , (ai) 

where x = r/e, with e denoted the softening length. This rule 
is optimal when particles are confined to a plane, because it 
yields in-plane accelerations that would result if the razor- 
thin mass distribution were displaced vertically by the soft- 
ening length. The forces can be thought of as approximating 
those from a disc of finite thickness since sof tening affects 
the dispersion relation for spiral waves (e.g. IVandervoortJ 
Il97d ; lEricksonlH975l : lRomedll992r ) in much the same way 
as does finite thickness. We therefore employ this rule when 
particle motion is confined to a plane. 

The disadvantage of the Plummer softening rule in 3D 
simulations is that it weakens forces on all scales and other 
rules that avoid this shortcoming have become popular. The 
precise short-range behaviour is of little importance for re- 
laxation, since inverse square-law forces imply scattering is 
dominated by the cumulative effect of long-range encoun- 
ters (e.g. BT08, p. 36). For our 3D simulations we adopt the 
somewhat clumsy cubic spline density kerne l used in th e 
original version of the tree code PKDGRAV l|Stadelll200il 'l. 
which has the form 

1 ( 4 - 6x 2 + 3x 3 < x < 1 
p( x ) = 7^3 \ (2 - z) 3 1 < x < 2 (A2) 

I otherwise, 

as suggested by iMonaghan fc Lattanziol (1985)0 The den- 
sity has continuous second derivatives, while the potential is 
given by a messy expression at short range but is, of course, 
that of a unit point mass when x > 2. 

A3 Codes 

For fully 3D simulat ions, we use t he hybrid grid method 
described elsewhere (|Sellwoodll2003l ). It solves for the field 
by convolution on a 3D cylindrical polar grid for the disc 
particles, with the softening rule (eq. IA2[I , while the acceler- 
ations of the halo particles are computed using method (3) 
of Section Al on a spherical grid. 

We also report a number of results using both polar and 
Cartesian 2D grids, where we use the Plummer softening 
law. 

In most experiments, we shift the centre of both grids 
to a new location every 16 tim e steps. The new centre is 
location of the particle centroid l|McGlvnnl l 19841'! . The esti- 
mate of the change in this location is determined by Newton- 
Raphson iteration, which is repeated until the shift at each 
iteration is less than 10~ 3 Rd- This process is unnecessary 
when any lop-sidedness in the mass distribution does not 
contribute to the accelerations and when Cartesian grids or 
tree codes are used. 

2 Their expression omits the square o f x in the first line, but this 
typo is corrected in iMonaehanl Jl992h . 
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In addit ion, we have used the tree code, PKDGRAV 
IStadell j200l[). which adop ts the softening kernel Ki recom- 
mended by iDehnenl {2001). We have also conducted a few 
tests with the numerical parameters of time step, opening 
angle, etc., and found results with this code also are inde- 
pendent of these choices to a similar level of tolerance. 



A4 Time steps 

When using grid methods, we adopt a 5-zone time stepping 
scheme in which the more slowly moving outer particles have 
time steps that increase by a factor 2 from zone to zone. 
All particles experience forces from all others at every step, 
but forces from par ticles in outer z ones are interpolated to 
intermediate times (|Sellwoodlll985l ). 



A5 Measurements of A and Q p 

We need to make quantitative comparisons of the bar evolu- 
tion in simulations as codes, numerical parameters, or ran- 
dom seeds are varied. In particular, we compare the evolu- 
tion of the overall amplitude and phase of a bar-like distor- 
tion in the disc. We measure this quantity by computing 



(A3) 



3=1 



where Nd is the total number of disc particles and 8j(t) is 
the azimuth of the jth disc particle at time t, reckoned from 
a fixed direction through the centre defined as the particle 
centroid at that time. Since the quantity A is complex, the 
bar amplitude is \A 2 \ = [%t(A 2 ) 2 + Q{A 2 ) 2 ] 1/2 and phase 
29b — arctan[9 ! (j4,2)/3i(A2)] ) with the factor 2 appearing 
in order to yield a phase that increases by 180° as the bi- 
symmetric pattern makes half a rotation. We measure A 2 at 
frequent intervals, generally every 0.1 dynamical times. The 
pattern speed of the bar is clearly the time derivative of the 
phase 

We make a smoothed estimate of the amplitude and 
pattern speed by fitting a steadily rotating wave to the com- 
plex A 2 values over a short time interval, and sliding the 
time interval forward to follow the evolution of both quan- 
tities. Our plots of amplitude and pattern speed are of the 
smoothed quantities. 



APPENDIX B: TESTS OF NUMERICAL 
PARAMETERS 

As always, we check the extent to which the behaviour de- 
pends upon all numerical parameters. We have been partic- 
ularly thorough in the case of our standard model where our 
results are so surprising. Since simulations with a rigid halo 
and the disc particles confined to a plane already show large 
variations (Fig. I15J1 , we begin by presenting checks of these 
inexpensive simulations. 



Bl Grid Geometry 

We have run these calculations on both a 2D polar grid and 
a 2D Cartesian grid in order to convince ourselves that our 





Figure Bl. Evolution of the bar in 16 runs with different random 
seeds for the disc particle coordinates. These simulations use a 
2D Cartesian grid: numerical parameters are as given in Table [2] 
except N x X N y = 256 2 , there is a single time step zone and the 
grid is not recentred. 

results were not being affected by our choice of grid geom- 
etry. The result for the Cartesian grid is shown in Fig. IB H 
which should be compared with that for the polar grid shown 
in the second row of Fig. 1151 for which the number of par- 
ticles and softening length were identical. Again the curves 
for separate runs have been shifted in time so that they all 
pass through amplitude 0.1 at the same instant, which is the 
mean of the set shown. 

While there are differences in detail between the two 
figures, the mean and spread in the amplitude evolution are 
quite similar. 

B2 Time Step 



Fig. IB2I shows that changing the time step also has little 
effect on the evolution. These tests are for two of the 3D 
models shown in Fig. [3] one in which the bar slowed at late 
times and one in which it did not. The value of the time 
step parameter is varied by a factor of 40 in the case that 
slowed strongly. Small differences in the evolution develop 
at late times because the system is chaotic but the devia- 
tions do not vary systematically with the step size. If the 
orbital angular frequency is Q c , a particle takes 2n / (Q, c At) 
steps for a circular orbit. The central value of fi c — 2 for our 
standard model, implying 250 steps per orbit for the most 
bound particles at our standard time-step, and ten times as 
man y for the shorte s t step used in Fig. IB2I In agreement 
with iDubinskT et al\ l|2009l), we therefore find no evidence 
to support the claim by Klypin et al\ (|2008T ) that these sim- 
ulations require > 2000 time steps per orbit period for the 
most tightly bound particles. 

We have also verified that the evolution is similarly in- 
sensitive to using a fixed time step for all particles, instead 
of the more efficient scheme of employing longer steps at 
larger distances from the centre. 

B3 Grid Resolution 

Fig. lB3l shows the effects of changing the size of the cylindri- 
cal polar grid used for the disc in the hybrid code, keeping 
the initial particle coordinates and all other numerical pa- 
rameters fixed. As in other tests, small differences in the evo- 
lution develop at late times but aside from the two coarsest 



© 2009 RAS, MNRAS 000.ITU2T1 



Stochasticity in N-body Discs 19 



0.4 I ' 1 ' i 1 0.5 



04 " 




— 0.006250 

0.012500 

0.025000 

0.050000 



200 400 600 800 ' 200 400 600 BOO 

time time 

Figure B2. Evolution of the bar in two sets of runs with the 
same random seeds for the disc particle coordinates. Numerical 
parameters are given in Table [T] except the time step is varied. 
Values adopted are colour coded as shown. Upper panels show the 
evolution from one set of initial coordinates, lower panels from a 
second set. 
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Figure B3. Evolution of the bar in 3D runs with the same ran- 
dom seeds for the disc particle coordinates using different grid 
sizes. Numerical parameters are given in Table [T] except the num- 
ber of grid cells used for the 3D polar grid is changed, as indicated 
by the line colours. 



grids, for which the late time evolution departs systemati- 
cally from the rest, the results are quite similar. We have 
also found that smaller differences result when we double 
or halve the vertical spacing of the grid. Our standard grid 
(Table [1]) is shown by the cyan line and seems adequate. 

In addition, we checked that the evolution is unaffected 
(to the same level of tolerance) by changing the number of 
active sectoral harmonics of the polar grid to m max = 4 
or m max = 16 from our standard value of m max = 8, or 
by changing the order of azimuthal expansion Z max and the 
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Figure B4. Evolution of the bar in 3D runs with the same sets 
of random seeds for the disc particles as in Fig. [4] but using dif- 
ferent softening lengths. The softening length in the upper and 
lower panels is respectively halved and doubled from our standard 
value. 

number of shells n r of the spherical polar grid from our 
standard values of Z max = 4 and n r = 300. 

B4 Softening 

Fig. IB4I shows the effects of changing the softening length 
for force convolution on the 3D polar grid used for the disc 
in the hybrid code, with the initial particle coordinates and 
other numerical p arameters unchanged from Fig. [3] As re- 
porte d elsewhere (|Sellwoodlll98ll . Il983l ; ISellwood fc MerrittJ 
1 19941 . etc.), the evolution of disc instabilities is more sen- 
sitive to this numerical parameter than perhaps any other. 
The growth rates of both bar-forming modes, and of bend- 
ing modes are quite sensitive to the sharpness of short-range 
forces. The effect of a longer softening length (lower panels) 
is both to increase the initial peak amplitude of the bar, 
because the second mode is more strongly suppressed by 
softening than is the dominant, and to make bending in- 
stabilities occur later and more violently0 The effects of a 
reduction in softening are less systematic, but the extra vir- 
ulence of swing-amplified shot noise is the probable cause 
of more marked upward fluctuations in the pattern speed 
evolution and there are fewer violent buckling events. 

Since it is desirable to use the largest value that does not 

3 The reason is as follows, adapted from iMerritt fe Sellwoodl 
||1994T) : If stars move at speed u in one-dimension over a rip- 
ple of wavenumber fc, then a condition for a growing bend is 
that ku < k z . Increasing softening reduces k z , which causes only 
smaller fc, or longer wavelength, bends to grow. 
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Figure B5. Evolution of the bar in runs to test the dependence 
on the number of halo particles. The upper panels used one tenth 
the number of unequal mass particles employed in Fig. [4] while 
the number employed in the lower panels was 10 million. Other 
numerical parameters are held fixed at the values given in Tablc[T] 



Figure CI. Evolution of the bar in the isochrone/8 disc, but in- 
stead of selecting particles deterministically as in Fig. 1131 we used 
a simple acceptance/rejection algorithm. Note the larger spread 
in the measured bar properties. 
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have a systematic influence on the outcome, these tests show 
that our standard value seems a reasonable compromise. 



Figure C2. Evolution of the bar in a noisy start isochrone disc 
in which the non-circular motions were set up crudely rather than 
selecting from a DF. The value of Q in the initial disc is similar 
to that of the initial models in Figs. Q2] Q3] felCTl 



B5 Number of Halo Particles 



Fig. IB5I shows two sets of runs with different numbers of 
unequal mass halo particles, in which the random seeds for 
the disc particles were changed. (We already reported the 
dependence of the behaviour on the number of disc parti- 
cles in Fig. 1151 ) Again, the behaviour in these tests, and in 
another set with 2.5M equal mass particles, is qualitatively 
similar to that shown in Fig [4] The ranges of final ampli- 
tudes and pattern speeds do not depend on the number of 
halo particles or whether the masses are all equal. There is 
a trend, in that the fraction of bars that do not experience 
strong friction seems to increase with increasing numbers 
of halo particles: it is 4/16 for N h = 2.5 x 10 5 , 7/16 for 
N h = 2.5 x 10 6 (Fig. g]| and 11/16 for N h = 10 7 . For the 
experiments with Nh = 2.5 x 10 6 equal mass particles, the 
non-slowing fraction is 3/16. 

We make use of this trend with the quality of the sim- 
ulations in the discussion of Sections 14.21 & [6] 



APPENDIX C: EFFECTS OF PARTICLE 
SELECTION FOR THE ISOCHRONE DISC 

Here we illustrate the advantages of careful particle selection 
for a simple disc model with well-defined global instabilities. 





Figure C3. Summary of results from in Figs, I13l f red). iMl f green). 
ICll fblue). &: IC2l ('cvan) in order to illustrate the ranges of scatter. 



The value of a quiet start was already illustrated by com- 
parison of Figs. [13] & [14] but particles were deterministically 
selected from the DF for both sets of simulations. 

Fig. IC1I shows the consequences of selecting particles 
by the commonly-used acceptance/rejection method. Even 
though these experiments still used quiet starts (replicas of 
each master particle spaced evenly around a ring), the re- 
sults are less well behaved: there is more scatter particularly 
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in the bar amplitude, with one or two significantly anoma- 
lous results. 

Fig. lC2l shows the results from experiments in which the 
set up procedure for the random speeds of the disc particles 
stemmed simply from the requirement that Q = 1.2 every- 
where, with the azimuthal dispersion and asymmetric drift 
determined by Jeans equations i n the epicycle approxima- 
tion, as suggested by iHernquistl |l993). Although this may 
be the most commonly used method, the outcome of such 
experiments shows the greatest degree of scatter. 

The effects of quiet and noisy starts, and other particle 
selection issues are summarized in Fig. IC3I Generally, ex- 
periments with noisy starts show considerably more scatter 
than do those with quiet starts, and deterministic selecting 
from a DF is superior to random sampling or not using a 
DF at all. 
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